#pragma once

#include <sfv/array1d.h>
#include <sfv/exception.h>
#include <sfv/detail/format_utils.h>

#include <map>
#include <tr1/random>
#include <time.h>

namespace sfv
{
	namespace gallery
	{
		template< typename IndexType>
		struct forest_fire_data_t 
		{
			typedef typename std::vector< sfv::array1d<IndexType> > Array;

			forest_fire_data_t( Array& inneis, Array& outneis, size_t no_of_nodes )
				: inneis(inneis), outneis(outneis), no_of_nodes(no_of_nodes){};

			Array& inneis;
			Array& outneis;
			size_t no_of_nodes;
		};

		/**
		* \function igraph_forest_fire_game
		* \brief Generates a network according to the \quote forest fire game \endquote
		* 
		* The forest fire model intends to reproduce the following network
		* characteristics, observed in real networks:
		* \ilist
		* \ili Heavy-tailed in-degree distribution. 
		* \ili Heavy-tailed out-degree distribution.
		* \ili Communities.
		* \ili Densification power-law. The network is densifying in time,
		*      according to a power-law rule.
		* \ili Shrinking diameter. The diameter of the network decreases in
		*      time.
		* \endilist
		* 
		* </para><para>
		* The network is generated in the following way. One vertex is added at
		* a time. This vertex connects to (cites) <code>ambs</code> vertices already
		* present in the network, chosen uniformly random. Now, for each cited
		* vertex <code>v</code> we do the following procedure:
		* \olist
		* \oli We generate two random number, <code>x</code> and <code>y</code>, that are
		*   geometrically distributed with means <code>p/(1-p)</code> and
		*   <code>rp(1-rp)</code>. (<code>p</code> is <code>fw_prob</code>, <code>r</code> is
		*   <code>bw_factor</code>.) The new vertex cites <code>x</code> outgoing neighbors
		*   and <code>y</code> incoming neighbors of <code>v</code>, from those which are
		*   not yet cited by the new vertex. If there are less than <code>x</code> or
		*   <code>y</code> such vertices available then we cite all of them.
		* \oli The same procedure is applied to all the newly cited
		*   vertices.
		* \endolist
		* </para><para>
		* See also:
		* Jure Leskovec, Jon Kleinberg and Christos Faloutsos. Graphs over time:
		* densification laws, shrinking diameters and possible explanations.
		* \emb KDD '05: Proceeding of the eleventh ACM SIGKDD international
		* conference on Knowledge discovery in data mining \eme, 177--187, 2005.
		* </para><para>
		* Note however, that the version of the model in the published paper is incorrect
		* in the sense that it cannot generate the kind of graphs the authors
		* claim. A corrected version is available from
		* http://www.cs.cmu.edu/~jure/pubs/powergrowth-tkdd.pdf, our
		* implementation is based on this.
		* 
		* \param graph Pointer to an uninitialized graph object.
		* \param nodes The number of vertices in the graph.
		* \param fw_prob The forward burning probability.
		* \param bw_factor The backward burning ratio. The backward burning
		probability is calculated as <code>bw.factor*fw.prob</code>.
		* \param pambs The number of ambassador vertices.
		* \param directed Whether to create a directed graph.
		* \return Error code.
		* 
		* Time complexity: TODO.
		*/

		template <typename MatrixType>
		int forest_fire_game(MatrixType& graph, int nodes,
							 double fw_prob, double bw_factor, int pambs) 
		{  
			typedef typename MatrixType::index_type IndexType;
			typedef typename MatrixType::value_type ValueType;

			sfv::array1d<long int> edges;
			long int no_of_nodes = nodes;

			sfv::array1d<long int> visited(no_of_nodes);
			std::deque<long int> neiq;
			long int ambs = pambs;
			double param_geom_out=1-fw_prob;
			double param_geom_in=1-fw_prob*bw_factor;

			if (fw_prob < 0) {
				throw sfv::invalid_input_exception("Forest fire model: 'fw_prob' should be between non-negative");
			}
			if (bw_factor < 0) {
				throw sfv::invalid_input_exception("Forest fire model: 'bw_factor' should be non-negative");
			}
			if (ambs < 0) {
				throw sfv::invalid_input_exception("Number of ambassadors ('ambs') should be non-negative");
			}

			if (fw_prob == 0 || ambs == 0) {
				throw sfv::invalid_input_exception("'fw_prob or ambs is zero, creating empty graph");
			}

			std::vector< sfv::array1d<IndexType> > inneis(no_of_nodes);
			std::vector< sfv::array1d<IndexType> > outneis(no_of_nodes);
			forest_fire_data_t<IndexType> data(inneis, outneis, no_of_nodes);

			srand((unsigned) time(NULL));

#define ADD_EDGE_TO(nei) \
	if (visited[(nei)] != actnode+1) {						\
	visited[(nei)] = actnode+1;								\
	neiq.push_back(nei);									\
	edges.push_back(actnode);								\
	edges.push_back(nei);									\
	outneis[actnode].push_back(nei);						\
	inneis[nei].push_back(actnode);							\
	}

			for (long int actnode=1; actnode < no_of_nodes; actnode++) 
			{
				/* We don't want to visit the current vertex */
				visited[actnode] = actnode+1;

				/* Choose ambassador(s) */
				for (long int i=0; i<ambs; i++) 
				{
					if( (actnode - 1) == 0 )
					{
						ADD_EDGE_TO(0);
					}
					else
					{
						long int a = rand() % (actnode-1);
						ADD_EDGE_TO(a);
					}
				}

				while (!neiq.empty()) 
				{
					long int actamb = neiq.front();
					neiq.pop_front();

					sfv::array1d<IndexType> & outv(outneis[actamb]);
					sfv::array1d<IndexType> & inv(inneis[actamb]);
					long int no_in = inv.size();
					long int no_out = outv.size();

					std::tr1::mt19937 mt;
					std::tr1::geometric_distribution<long int, double> geometric_out(param_geom_out); 
					long int neis_out = geometric_out(mt);
					std::tr1::geometric_distribution<long int, double> geometric_in(param_geom_in); 
					long int neis_in = geometric_in(mt);

					/* outgoing neighbors */
					if (neis_out >= no_out) {
						for (long int i=0; i<no_out; i++) {
							long int nei = outv[i];
							ADD_EDGE_TO(nei);
						}
					} else {
						long int oleft=no_out;
						for (long int i=0; i<neis_out && oleft > 0; ) {
							long int which = rand() % (oleft-1);
							long int nei = outv[which];
							outv[which] = outv[oleft-1];
							outv[oleft-1] = nei;
							if (visited[nei] != actnode+1) {
								ADD_EDGE_TO(nei);
								i++;
							}
							oleft--;
						}
					}
					/* incoming neighbors */
					if (neis_in >= no_in) {
						for (long int i=0; i<no_in; i++) {
							long int nei = inv[i];
							ADD_EDGE_TO(nei);
						}
					} else {
						long int ileft=no_in;
						for (long int i=0; i<neis_in && ileft > 0; ) 
						{
							long int which;
							
							if( (ileft - 1) == 0 )
							{
								which = 0;
							}
							else
							{
								which = rand() % (ileft-1);
							}
							long int nei = inv[which];
							inv[which] = inv[ileft-1];
							inv[ileft-1] = which;
							if (visited[nei] != actnode+1) {
								ADD_EDGE_TO(nei);
								i++;
							}
							ileft--;
						}
					}

				} /* while neiq not empty */

			} /* actnode < no_of_nodes */

			graph.resize( no_of_nodes, no_of_nodes, edges.size() );
			sfv::array1d<IndexType> row_indices(edges.size());

			std::multimap<int,int> edge_map;

			for(size_t i = 0; i < edges.size()-1; i+=2 )
			{
				edge_map.insert(std::pair<int, int>(edges[i], edges[i+1]));
				edge_map.insert(std::pair<int, int>(edges[i+1], edges[i]));
			}

			size_t index = 0;
			for (std::multimap<int,int>::iterator it = edge_map.begin(); it != edge_map.end(); ++it)
			{
				row_indices[index] = it->first;
				graph.column_indices[index++] = it->second;
			}

			sfv::detail::indices_to_offsets(row_indices, graph.row_offsets);

			return 0;
		}


	} // end namespace gallery
} // end namespace sfv
